"""
Cahn-Hilliard PDE implementation

The Cahn-Hilliard equation describes phase separation and spinodal decomposition:
    ∂c/∂t = ∇²(c³ - c - γ∇²c)

This is a fourth-order PDE that conserves the order parameter c.
"""

using ..Fields
using ..Grids

"""
    CahnHilliardPDE <: AbstractPDE

Cahn-Hilliard equation for phase separation and spinodal decomposition.

The mathematical definition is:
    ∂c/∂t = ∇²(c³ - c - γ∇²c)

where c is the order parameter (concentration field) and γ controls the interface width.

# Parameters
- `interface_width`: Controls interface width γ (default: 1.0)
- `mobility`: Overall mobility parameter (default: 1.0)
- `bc_c`: Boundary conditions for the order parameter c
- `bc_mu`: Boundary conditions for the chemical potential μ

# Properties
- Fourth-order PDE requiring special handling
- Conserves total mass: ∫c dV = constant
- Energy is non-increasing: dE/dt ≤ 0
- Exhibits spinodal decomposition for |c| < 1
"""
struct CahnHilliardPDE <: AbstractPDE
    interface_width::Float64
    mobility::Float64
    bc_c::String
    bc_mu::String
    
    function CahnHilliardPDE(; interface_width=1.0, mobility=1.0, 
                            bc_c="auto_periodic_neumann", bc_mu="auto_periodic_neumann")
        new(interface_width, mobility, bc_c, bc_mu)
    end
end

"""
    evolution_rate(pde::CahnHilliardPDE, state::ScalarField, t::Float64)

Compute the evolution rate for the Cahn-Hilliard equation.

The equation is: ∂c/∂t = ∇²μ where μ = c³ - c - γ∇²c
This is implemented as two steps:
1. Compute chemical potential: μ = c³ - c - γ∇²c
2. Apply mobility: ∂c/∂t = M∇²μ
"""
function evolution_rate(pde::CahnHilliardPDE, state::ScalarField, t::Float64)
    c = data(state)
    grid_obj = grid(state)
    
    # Step 1: Compute Laplacian of c
    c_laplace = laplacian(state)  # This returns a ScalarField
    
    # Step 2: Compute chemical potential μ = c³ - c - γ∇²c
    # Create cubic nonlinearity field
    c_cubic_data = c.^3
    c_cubic = ScalarField(c_cubic_data, grid_obj)
    
    # Create linear term field  
    c_linear = ScalarField(-c, grid_obj)
    
    # Chemical potential: μ = c³ - c - γ∇²c
    mu = c_cubic + c_linear - pde.interface_width * c_laplace
    
    # Step 3: Apply mobility and compute ∇²μ
    mu_laplace = laplacian(mu)
    
    # Return M∇²μ
    return pde.mobility * mu_laplace
end

"""
    free_energy(pde::CahnHilliardPDE, state::ScalarField)

Compute the free energy functional for the Cahn-Hilliard system:
    E = ∫[f(c) + (γ/2)|∇c|²] dV
where f(c) = (c² - 1)²/4
"""
function free_energy(pde::CahnHilliardPDE, state::ScalarField)
    c = data(state)
    
    # Bulk free energy density: f(c) = (c² - 1)²/4
    bulk_energy = sum((c.^2 .- 1).^2) / 4
    
    # Gradient energy density: (γ/2)|∇c|²
    grad_c = gradient(state)
    grad_magnitude_sq = sum(norm.(grad_c)^2 for grad_c in grad_c.data)
    gradient_energy = pde.interface_width * grad_magnitude_sq / 2
    
    return bulk_energy + gradient_energy
end

"""
    is_mass_conserved(pde::CahnHilliardPDE)

Check if the PDE conserves mass (total integral of c).
Returns true for Cahn-Hilliard equation.
"""
is_mass_conserved(pde::CahnHilliardPDE) = true

"""
    get_characteristic_length(pde::CahnHilliardPDE)

Get the characteristic interface width for the Cahn-Hilliard equation.
"""
get_characteristic_length(pde::CahnHilliardPDE) = sqrt(pde.interface_width)

"""
    get_time_scale(pde::CahnHilliardPDE, grid_spacing)

Estimate the characteristic time scale for the Cahn-Hilliard dynamics.
For stability, dt should be much smaller than this.
"""
function get_time_scale(pde::CahnHilliardPDE, grid_spacing)
    # Time scale is roughly L⁴/(M·γ) where L is grid spacing
    return grid_spacing^4 / (pde.mobility * pde.interface_width)
end

"""
    is_stiff(pde::CahnHilliardPDE)

The Cahn-Hilliard equation is generally stiff due to fourth-order derivatives.
"""
is_stiff(pde::CahnHilliardPDE) = true

"""
    requires_high_accuracy(pde::CahnHilliardPDE)

The Cahn-Hilliard equation typically requires high accuracy to capture
interface dynamics correctly.
"""
requires_high_accuracy(pde::CahnHilliardPDE) = true

"""
    expression(pde::CahnHilliardPDE)

Return a string representation of the PDE.
"""
function expression(pde::CahnHilliardPDE)
    γ_str = pde.interface_width == 1.0 ? "" : "$(pde.interface_width)"
    M_str = pde.mobility == 1.0 ? "" : "$(pde.mobility) "
    if isempty(γ_str)
        return "$(M_str)∇²(c³ - c - ∇²c)"
    else
        return "$(M_str)∇²(c³ - c - $(γ_str)∇²c)"
    end
end

# Export the PDE type
export CahnHilliardPDE, free_energy, is_mass_conserved, get_characteristic_length, get_time_scale